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ABSTRACT 

Measurement of the angular power spectrum of the cosmic microwave background 
is most often based on a spherical harmonic analysis of the observed temperature 
anisotropics. Even if all-sky maps are obtained, however, it is likely that the region 
around the Galactic plane will have to be removed due to its strong microwave emis- 
sions. The spherical harmonics are not orthogonal on the cut sky, but an orthonormal 
basis set can be constructed from a linear combination of the original functions. Pre- 
vious implementations of this technique, based on Gram-Schmidt orthogonalisation, 
were limited to maximum Legendre multipoles of /,nax ^ 50 as they required all the 
modes have appreciable support on the cut sky, whereas for large /max the fraction 
of modes supported is equal to the fractional area of the region retained. This prob- 
lem is solved by using a singular value decomposition to remove the poorly-supported 
basis functions, although the treatment of the non-cosmological monopole and dipole 
modes necessarily becomes more complicated. A further difficulty is posed by compu- 
tational limitations - orthogonalisation for a general cut requires 0{lf^^^^) operations 
and 0(/max) storage and so is impractical for l^ax ^ 200 at present. These problems 
are circumvented for the special case of constant (Galactic) latitude cuts, for which the 
storage requirements scale as 0(/max) ^^'^ the operations count scales as 0{lfj^^^). Less 
clear, however, is the stage of the data analysis at which the cut is best applied. As 
convolution is ill-defined on the incomplete sphere, beam-deconvolution should not be 
performed after the cut, and, if all-sky component separation is as successful as simu- 
lations indicate, the Galactic plane should probably be removed immediately prior to 
power spectrum estimation. 

Key words: cosmic microwave background - methods: analytical - methods: numer- 
ical. 



1 INTRODUCTION 

Since the first measurements of the temperature anisotropy 
of the cosmic microwave background (CMB) by the Cos- 
mic Background Explorer (COBE) satellite (Smoot et al. 
1992), a number of sophisticated experiments have been un- 
dertaken to measure the fluctuations at higher resolutions 
and sensitivities (e.g. Scott et al. 1996; Tanaka et al. 1996; 
Netterfield et al. 1997; de Oliveira-Costa et al. 1998; Coble 
et al. 1999; de Bernardis et al. 2000; Wilson et al. 2000; 
Padin et al. 2001; Halverson et al. 2001; Lee et al. 2001; 
Netterfield et al. 2001). The primary result of these experi- 
ments has been the measurement of the angular power spec- 
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trum of the CMB to Legendre multipoles of up to / ~ 1000, 
which places strong constraints on a number of cosmologi- 
cal parameters (Lineweaver 1998; Efstathiou et al. 1999; de 
Bernardis et al. 2000; Netterfield et al. 2001; Wang, Tegmark 
& Zaldarriaga 2001 and references therein). In the future the 
Microwave Anisotropy Probe {MAP; e.g. Jarosik et al. 1998) 
and the Planck satellite (e.g. Bersanelli et al. 1996) will pro- 
duce maps of the microwave sky with resolutions of between 
5 and 30 arcmin at a number of frequencies. Such extraor- 
dinary data-sets, consisting of millions of independent mea- 
surements, will clearly require novel analysis techniques. 

One of the many difficulties is the treatment of the non- 
cosmological contributions to the observed microwave sky. 
Dust, synchrotron and free- free emission from the Galaxy 
(e.g. Haslam et al. 1982; Schlegel, Finkbinder & Davies 
1998); radio galaxies and other extra- Galactic 'point' sources 
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(e.g. TofFolatti et al. 1998); and the Sunyaev-Zel'dovich 
(1970) effect caused by galaxy clusters (e.g. Birkinshaw 
1999) all obscure the CMB at some level (see Hu, Sugiyama 
& Silk 1997 or Barreiro 2000 for more complete reviews), al- 
though these components have quite distinct spectral prop- 
erties and so can be separated using multi-frequency ob- 
servations (e.g. Bennett et al. 1992; Tegmark & Efstathiou 
1996; Hobson et al. 1998; Bouchet & Gispert 1999; Jones, 
Hobson & Lasenby 1999; Baccigalupi et al. 2000). However 
these techniques are not likely to be able to extract the 
Galactic emissions completely (Stolyarov et al. 2001), leav- 
ing the removal of the Galactic plane as the only option. The 
Galaxy contributes relatively little at high latitudes (e.g. 
Haslam et al. 1982; Schlegel et al. 1998) so this is an ac- 
ceptable, if not optimal, solution. For instance, Gorski et al. 
(1994) removed the band within 20 deg of the Galactic plane 
to estimate the power spectrum of the two-year COBE Dif- 
ferential Microwave Radiometer sky maps, and similar cuts 
have been proposed by both the MAP and Planck collab- 
orations. An essentially equivalent problem is posed if the 
survey's sky coverage is incomplete, although there is less 
choice about the geometry of the cut in this case. 

A number of aspects of the analysis become more diffi- 
cult on an incomplete sphere, one of the most obvious rea- 
sons being that the spherical harmonics are no longer an 
orthonormal basis set. The most successful component sep- 
aration techniques to date (e.g. Hobson et al. 1998; Bouchet 
& Gispert 1999) rely on a mode-by-mode analysis which ex- 
plicitly utilises the orthogonality of the spherical harmonics, 
although it may be preferable to remove the Galactic plane 
only when estimating the CMB power spectrum. Unbiased 
power spectrum estimation using the spherical harmonics is 
possible on the cut sky (Wandelt, Gorski & Hivon 2001), but 
the covariance structure of the resulting psuedo-harmonics 
is far from ideal, so analysis using an orthonormal basis set is 
preferable. In particular the noise covariance matrix remains 
diagonal in the case of spatially uniform (and uncorrelated) 
noise. 

It is possible to construct an orthonormal basis set from 
linear combinations of the spherical harmonics, and an ele- 
gant implementation of this, based on Cholesky decomposi- 
tion of the coupling matrix of the spherical harmonics on the 
cut sky, was described by Gorski (1994). However the cou- 
phng matrix becomes ill-conditioned for /max ^ 50, and so 
this method cannot be used to perform cut-sky orthogonal- 
isation for either the MAP experiment (with imax — 1500) 
or the Planck mission (with /max — 2500). 

A general formalism for orthogonalisation of the spheri- 
cal harmonics is presented in Section ^, although implemen- 
tation to high iniax is only possible at present in the special 
case of a constant latitude cut. The relationship between the 
various harmonic coefficients is discussed in Section ^, and 
the extension of these results to CMB analysis techniques 
(specifically component separation and power spectrum es- 
timation) is covered in Section ^. The results are summarised 
and future possibilities are discussed in Section H Finally, 
the chosen conventions for the spherical harmonics are de- 
fined in Appendix^; formulae for integrals of the products of 
Legendre functions are given in Appendix ^ and the treat- 
ment of the non-cosmological monopole and dipole modes 
are discussed in Appendix ^. 



2 ORTHOGONALISATION OF SCALAR BASIS 
FUNCTIONS 

The physics of the CMB is most naturally expressed in 
Fourier space, and it is standard practice to represent sky 
maps by their harmonic coefficients. The basis functions cho- 
sen here are the real spherical harmonics, Yi^mir) (as defined 
in Appendix which form an orthonormal basis on the 
complete sphere, S. In general I > and —l<m<l, al- 
though in practice a finite imax must be used, which im- 
plies a band-limit. It is convenient to combine the two 
indices, allowing the basis set to be expressed as a vec- 
tor, Y{r) = [Vi(r),y2(f-),...,y>„_(f)]'^, where w = 
(/max + 1)^. There are several reasonable choices for the in- 
dexing, i{l,m), most notably grouping coefficients in / or 
m, as defined in Appendix Ul. Grouping in I is most nat- 
ural for power spectrum estimation, but grouping in m is 
more efficient com putationally in cases of azimuthal sym- 
metry (Section 2.3). 

The spherical harmonics are not orthogonal on the in- 
complete sphere S' , as ca n be seen from the structure of their 
coupling matrix (Section 2.1). A decomposition of the cou- 
pling matrix can be used to construct an orthonormal basis 
set (Section ^.2| ), but implementation to high-resolution is 
currently possible only in the special case of constant lati- 
tude cuts (Section) 



2.1 The coupling matrix 

The coupling matrix of a set of functions encodes their or- 
thogonality and normalisation properties over a given range. 
In the case of the spherical harmonics on the incomplete 
sphere it is given by 



C 



Y{r)Y'^{r)AQ.. 



(1) 



\i S' = S then the harmonics are orthonormal and C = I; 
otherwise the off-diagonal elements are non-zero, indicating 
that the basis functions are non- orthogonal. 

An alternative formulation is to introduce a window 
function, w{r), so that 

C = I u?{r)Y{r)Y'^{r)An. (2) 

JS 

In some ways this approach is more flexible, as w{r) can ei- 
ther be a smoothly-varying apodizing function (cf. Tegmark 
1997) or take the form 



ws'(r) = 



1, if f G S", 
0, if f£S-S', 



(3) 



mimicing the effect of the sharp cut deflned above. However 
this definition of the window function can lead to inconsis- 
tencies if a band-limited analysis is carried out, as ws'ir) 
cannot be properly represented by a finite analysis (see Sec- 
tion |j.2| ). It is for this reason that the first formalism is used 
here, although most of the subsequent results can also be 
derived using window functions. 

For a pixel-based analysis, the coupling matrix can be 
defined by replacing the integral in equation (jl|) by a sum 
over points on the sphere (i.e. pixel centres), fp, where 1 < 
p < Np and A'^p is the number of pixels. In this case 
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Afp 

C= Y{r^)Y^{rp)n^, (4) 

p=l;rpGS' 

where Q.p is the area of the pth pixel and there is no pixel- 
smoothing (cf. Gorski 1994). In the limit A^p oo equa- 
tions (0) and become equivalent and pixelisation issues 
become irrelevant. If the points are uniformly distributed 
over the sphere C should be close to the identity, the small 
discrepancies merely reflecting the approximation of the in- 
tegral as a sum; otherwise C reflects the spatial distribution 
of the points much as in the continuum case, but there is 
freedom to represent apodizing filters as well as discrete cuts. 

Assuming fig/ > 0, the coupling matrix is formally sym- 
metric, positive definite and invertible, irrespective of which 
of the above definitions is used. However C rapidly becomes 
numerically singular: e.g. if Zmax = 50, the condition num- 
bei^of C is ~ 5x 10® for a constant latitude cut of &cut = ±20 
deg. This can be further understood in terms of the eigen- 
structure of the coupling matrix. 



corresponding to the null-space of C) with the Ai = 1 modes 
(i.e. those in the range of C), and so the number of supported 
modes is determined by a combination of the band-limit and 
the cut. 

The number of orthonormal basis functions (i.e. the 
rank of C) is proportional to the area of the sphere retained, 
0.$' ■ Hence it is possible to define only 

i' ~ - rs) 

'-max — f-^ tniax — . ^max \<J } 

\ls 47r 

orthonormal functions on the cut sphere for a given (large) 
band-limit. The relative reduction in the basis set is the 
same as would occur in the equivalent pixel analysis: the 
number of pixels retained is also given by A'^^ ~ Q,si /^s Np. 
For low /max thcse arguments do not hold, and it is possible 
to create a basis set with more than flgi /Qs imax elements. 
Moreover, all these functions are required to ensure that the 
cut-sphere basis set is complete (as well as orthonormal) in 
the case of a low band-limit. 



2.1.1 Eigenstructure 

The coupling matrix has imax = (/max -I- 1)^ eigenvectors, Vi, 
and eigenvalues, Xi, which satisfy 



Cvi = XiVi. 



(5) 



Premultiplying by Y'^{r) and expanding out the implicit 
summations gives 

^max *max 

£ Yk{r)iv,)k 2 Yj(r)YAr) dil' = XMr), (6) 
fc=i i=i 

where Vi{r) = Y'^{r)vi is the ith eigenfunction of the cou- 
pling matrix. The completeness of the spherical harmonics 
in the /max ^ oo limit implies that Y"^ {r)Y {r') = 5{r ~r' 
[where S{x) is the Dirac delta function], so that equation i 
reduces to 



Vi{r')S{r — r) dQ' = XiVi{r). 



(7) 



For a given i this must be true at all f, which implies that 
either: Vi{r) = in the cut region, S — S' , in which case 
Xi — 1; or Vi(r) = in S' , in which case Xi — 0. In other 
words these eigenfunctions are completely localised in either 
the cut sphere or the removed region. This bimodality is only 
strictly true in the /max — * oo limit, but, as shown in Fig. nl 
is a good approximation for /max ^ 500. 

As the coupling matrix is symmetric, those eigenvectors 
with different eigenvalues are orthogonal, and those with 
the same eigenvalues can be made orthogonal by a rotation 
in the subspace defined by the eigenvalue in question (e.g. 
Arfken 1985). Thus the eigenfunctions with Xi = 1 represent 
an orthogonal basis set on S' , whereas those with Ai = 
have no support in this region and so cannot be orthogonal 
(or normalised) on the cut sky. The freedom in choice of 
basis does not extend to mixing the Xi = modes (i.e. those 



t The condition number of a matrix is the (absolute value of) 
the ratio of its greatest and smallest eigenvalues; it is large for 
ill-conditioned matrices, and infinite for singular matrices. 



2.2 Construction of an orthonormal basis 

The construction of an orthonormal basis set from a set of 
linearly-independent functions is a well-established mathe- 
matical technique, and a number of orthogonalisation meth- 
ods are possible. The most basic is Gram-Schmit orthogo- 
nalisation (e.g. Arfken 1985), in which the new basis func- 
tions are built-up sequentially, but this algorithm is numeri- 
cally unstable. The modified Gram-Schmidt algorithm (e.g. 
Golub & van Loan 1996) is stable, but it is generally prefer- 
able to use matrix techniques to create all the new basis 
functions simultaneously. 

Starting with the spherical harmonics, Y{r), the task 
is to find a set of functions Y'{r) which are orthonormal on 
the incomplete sphere. In terms of a conversion matrix, B, 
the two sets of functions are related by 

Y'{r) = BY{r). (9) 

Note that B has dimensions imax x imax, where imax ~ 
(/max + 1)^ and imax 1^ *max is determined by the band-limit 
and the cut, as described in Section It is also impor- 
tant to note that the indexing of the Yi (f) is qualitatively 
different from the Yi{r). The latter are really two-index ob- 
jects, with their characteristic scale given by ~ tt// and m 
relating to 'orientation'. However the new basis functions in- 
clude contributions from spherical harmonics with different 
/-values, and thus do not have a well-defined angular scale. 
Hence their single index contains no physical information, 
and the ordering or grouping of the new basis functions is 
arbitrary. 

From equation (|^), the coupling matrix of these new 
functions is 

c' = / Y'(f)y''^(f)do 

Js' 

= BCB'^. (10) 
Hence any conversion matrix which satisfies^ 



t Here I is the imax ^ * 
(potentially larger) i 



max identity matrix, as distinct from the 
max X imax identity matrix, I. 
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Figure 1. The distribution of eigenvalues of C, shown for symmetric constant latitude cuts of fecut = 10 deg in (a) and fecut = 20 deg in 
(b). In both panels the eigenvalues are sorted into decreasing order, and distributions are shown for Zmax = 10 (solid line), /max = 100 
(dashed line) and /max = 1000 (dotted line). The limiting cases of /max — ► oo [in which S73//(47r) of the eigenvalues are unity and the 
remaineder zero] are almost indistinguishable from the /max = 1000 distributions. 



bcb"^ = I 



(11) 



2.2.1 Low-resolution analysis 



yields basis functions which are orthonormal on the cut 
sphere, and the task of orthogonalisation is reduced to find- 
ing a solution for B given C. Whilst such a solution does 
not exist for arbitrary C, in all cases of practical interest 
a suitable conversion matrix can be constructed from the 
coupling matrix. One possible method is direct calculation 
of the eigenstructure of C, which yields a conversion matrix 
with elements given by Bii i = (vii)i\~,^^'^ , where the Vii 
are the eigenvectors of C and the A;/ its positive eigenval- 
ues. However it is advantageous to include the symmetry of 
the coupling matrix explicitly, which leads to a factorization 
of the form 



AA' 



(12) 



termined by the decomposition method. Combining equa- 
tions ( prj ) and (p^, the task of orthogonalisation is reduced 
to finding B such that 



BA = 0, 



(13) 



orthogonal matrix (i.e. 00 



where is an ij^^x 
= I). 

Whilst equations (|l^) and (^^ are general expressions 
which must be satisfied by the conversion matrix, they do 
not define a definite algorithm for the orthogonalisation. In 
practice it is simplest to choose = 1, leading to the re- 
quirement that 



BA = I. 



(14) 



However the optimal choice of decomposition method used 
to generate A depends on whether the coupling matrix is 
(numerically) invertible, and hence on the band-limit of the 
analysis. 



If 'max Si 50 and most of the sphere is retained (i.e. fig' ^ 
f7s/2 = 2ii) then the coupling matrix is numerically in- 
vertible and can be treated as positive definite in practice. 
Consequently A [defined in equation ([L^] is also invertible, 
and B — A~^, so that, from equation^l), the orthonormal 
basis set is given by 



A-'Y(f) 



(15) 



The form of A depends on the factorization method; of the 
wide variety available (e.g. Golub & van Loan 1996), the two 
most useful here are singular value decomposition (SVD) 
and Cholesky decomposition. 

The SVD of the covariance matrix is defined in terms 
of equation (|l|) by A = VW^/^ (i.e. C = VVVV"^), where V is 
orthogonal and W is diagonal^. The diagonal elements of W 
are the eigenvalues of C and, as their ordering is arbitrary, 
W can be defined such that Wi^i < Wi+i^i+i, provided the 
columns of V are permuted in the same way. The columns 
of V, in turn, are the eigenvectors of C, and V is an orthogo- 
nal matrix (i.e. = V^)- Hence the conversion matrix is 
given by B = A"'^ — W~^''^V^, which is trivially computed 
once the SVD has been performed. Note that this approach 
is effectively the same as the direct calculation of the eigen- 



structure of C mentioned above in Section 2.2 

Whilst SVD is a powerful technique, it is computation- 
ally expensive - a Cholesky decomposition is approximately 
ten times faster, although it can only be performed on sym- 
metric matrices which are numerically positive definite. The 
Cholesky decomposition of the covariance matrix takes the 
form C = LL^ (i.e. A = L), where L is lower triangular. Hence 



§ If M is diagonal then the notation M^^''^ is used here to de- 
note the matrix defined by (M^^/'^)ij = SijM^^'''^ , where Sij 

is the Kronecker delta function. Thus M^/'^ only exists if the di- 
agonal elements of M are non-negative and M~^^^ only exists if 
the diagonal elements of M are strictly positive. 
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Figure 2. Orthonormal cut sphere basis functions, Yl{r), as described in Section 2.2.2. In all cases imax = 100, W^in = 0.01, and a 
MoUweide projection is used. The colour-map varies from black (large negative values), through grey (zero), to white (large positive 
values), in each case being scaled to cover the dynamic range of the relevant basis function. For those in the left column [(a), (c), (e) and 
(g)] the cut (shown by the dashed lines) is symmetric, with f)cut = ±20 deg; for those in the right column [(b), (d), (f) and (h)] the cut 
(again shown by the dashed lines) is asymmetric, the region between fei = deg and 62 = —45 deg having been removed. Within each 
column the indexing of the basis functions is arbitrary, but they are displayed so that their fractional support in the removed region 
increases from (a) to (g) and (b) to (h), respectively. 
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the conversion matrix, B = can be computed quickly 
from the initial factorization in this case as well. The trian- 
gular structure of the conversion matrix also ensures that 
the new basis functions are the same as those formed by a 
numerically-stable Gram-Schmidt orthogonalisation (Gorski 
1994). 

Despite the fact that the SVD and the Cholesky de- 
composition result in quite different sets of basis functions, 
there is no reason to prefer one over the other in general. In 
the case of CMB analysis, however, the triangular structure 
of A and B as generated by the Cholesky decomposition is 
preferable as it ensures that the non-cosmological monopole 
and dipole modes are kept separate from the I > 2 modes, 
assuming /-ordering is used (Gorski 1994). If the SVD route 
is taken (or another indexing scheme used) the separation of 
the I — and I = 1 modes can be ensured usingthe partial 
Householder transform described in Appendix O Nonethe- 
less, if the coupling matrix is sufficiently non-singular, a 
Cholesky decomposition should be used to create the or- 
thonormal basis set, due to both its computational efficiency 
and the simplicity with which the non-cosmological modes 
are handled. 



2.2.2 High-resolution analysis 

If /max ^ 50 the couphng matrix is numerically singular, and 
thus A [defined in equation (^2|)] is non-invertible. Cholesky 
decomposition of C is thus impractical and, whilst an SVD 
is possible, the conversion matrix as defined in Section 



2.2.1 



cannot be computed, as the smallest elements of W (i.e. the 
smallest eigenvalues of C) are so close to zero. This implies 
that the corresponding columns of V do not contribute to 
the reconstruction of C and can be ignored. Hence it is possi- 
ble to perform an approximate SVD of the coupling matrix, 
defined by C ~ VVVv"^ (i.e. A = VW^''^), where W is an 
*max X iniax diagonal matrix containing the largest elements 
of W and V is an i^^ax ^ *max matrix consisting of the corre- 
sponding columns of V. The value of imax is determined by 
the choice of Wmin used to truncate W, but the bimodality 
of the eigenvalue distribution means that any value between 
~ 10~^ and ~ 0.1 is acceptable. The resultant conversion 

~— 1/2~T ~ 

matrix is B = W V (satisfying BA ~ I) and the i'^^ax 
new basis functions are given by 



Y'{r) 



w '''^v^y(f). 



(16) 



These basis functions represent an orthonormal basis set 
on the incomplete sphere, but they are not formally complete 
to the nominal band-hmit due to the slightly approximate 
nature of the reduced SVD. The decomposition becomes ex- 
act in the limit /max ^ oo as W ^ I (from the eigenvalue 



arguments described in Section 2.1.1) and the reduced SVD 

— x ~ ' ' ~ T 

becomes C VV (i.e. A — > V and B ^ V ) in this limit. 

Note also that these basis functions are still orthogonal on 

the full sphere, although they are no longer normalised to 

unity; this is another potential advantage of the SVD-based 

method. 

Several examples of orthonormal cut-sphere basis func- 
tions are shown in Fig. |^, for both symmetric and asymmet- 
ric constant latitude cuts. The link between these functions 
and the spherical harmonics is apparent - they have the 
same cellular structure, but, for the most part, are com- 



bined in such a way that their support is localised to S' . 
However the functions shown in Fig. ^ (g) and (h) are from 
the small subset with intermediate values of Wi^i and, as 
such, have considerable support in the removed region. 

In the case of CMB analysis, one short-coming of this 
approach is that the non-cosmological monopole and dipole 
modes are not distinguished from the higher moments (see 
Section 3.2), although this separation can be achieved post 
facto by usinga partial Householder transform, as described 
in Appendix O. In doing this some of the useful properties 
of the SVD are lost, but this operation need only be per- 
formed as the last step in generating the orthonormal basis 
set, by which stage all the computationally intensive matrix 
operations have already been performed. 

Despite the inconvenience caused by the mixing of the 
monopole and dipole modes, SVD is clearly the most fiexible 
and general method of orthogonalisation. In part, this stems 
from the fact that it can be applied to the coupling matrix 
without any prior knowledge of its singular properties. How- 
ever, in the high-/max limit, the number of supported basis 
functions is determined by the area of the cut and given, 
to a good approximation, by equation (^. Hence faster, if 
less powerful, techniques can be used to orthogonalise the 
spherical harmonics for /max <; 1000 as the coupling matrix 
is guaranteed to have at least (/max + l)'^ns'/(47r) positive 
eigenvalues. An example of this idea would be to modify the 
pivoting algorithm in the psuedo-ChoIesky decomposition 
described in Section 4.2.9 of Golub & van Loan (1996) so 
that the decomposition halts when the predetermined num- 
ber of basis functions have been generated, as opposed to 
using the less robust threshold based on the values of the 
diagonal elements of C. 



2.3 Constant latitude cuts 

In principle the method presented above is a complete so- 
lution to the problem of constructing orthonormal bases on 
the cut sphere, but the coupling matrix requires 0(/^ax) 
storage, limiting a general implementation to /max — 200 
on most current computers. Furthermore, the SVD of an 
n X n matrix requires O(n^) operations, and so the orthog- 
onalisation operation count scales as ©(/max)- Similar dif- 
ficulties are encountered in merely evaluating C, regardless 
of whether numerical integration or recursive techniques are 
used (e.g. Hivon et al. 2001). 

Fortunately all these difficulties are significantly re- 
duced in the case of a constant latitude cut (cf. Oh, Spergel 
& Hinshaw 1999; Wandelt et al. 2001), defined by ignoring 
all 9 for which 6\ < 9 < 62. This could be the symmet- 
ric removal of the Galactic plane (i.e. 6\ = 7r/2 — ?)cut and 
62 = 7r/2 -I- /)cut, where fecut is the latitude of the cut) or the 
absence of data round one pole (i.e. Si = and 62 = Scut). 
The formalism derived below can also be trivially extended 
to include multiple cuts, as would be required for a CMB 
experiment which did not observe either ecliptic pole. 

Explicitly including the constant latitude cut in equa- 
tion (^), the elements of the coupling matrix are given by 



C, 



1(1,171) ,m' ) 



1 



(17) 
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cos(f^2) 



+ 



cos{^?l ) 



where Sm{<t>) is defined in equation (A2), and the \i,m(x) 
are n orma hsed associated Legendre functions, given in equa- 
tion (|A3| ). From Appendix |A|, the first integral in equa- 
tion (O) reduces to 27r5m and so 



^ i(l ,vn) ,m' ) ^m,m' 



(18) 



/■cos(ei) 

<5i,i' - 27r / A,,|„|(a::)A,/, 1^1(2::) da; 

Jcos{92) 

The remaining integral can be evaluated using a combina- 
tion of analytical formula and recursion relations, as de- 
scribed in Appendix ^. 

The most important aspect of equation ( [l^ ) is that the 
coupling matrix C is extremely sparse (only one element in 
~ ^max is non-zero) a nd, if stored using the indexing scheme 
defined in equation (A.12) (i.e. grouped into sub- matrices 
of fixed m), is block diagonal. C can thus be stored in the 
form of 2/max + 1 sub-matrices, the mth of which has (/max + 
1 — Im])^ elements, and the storage requirements thus scale 
as 0(ij5iax) rather than 0(/J^ax)- Whilst it is convenient to 
store all the blocks simultaneously, there is no need to do 
so, which can further reduce the storage requirements to 
0(/max)- It is also clear from equation (^) that only the 
m > terms need be treated explicitly and that / and I' 
are interchangeable, decreasing the storage requirements by 
a further factor of four. Finally, in the case of a symmetric 
cut (i.e. O2 — n — 6\) the parity of \i,m(x) is such that all 
terms for which l + l' is odd vanish, resulting in an additional 
halving of the memory requirements. 

The orthogonalisation can be performed by decom- 
posing each sub-matrix separately, reducing the operation 
count from O(Z^ax) to O(Zmax)- The removal of the poorly- 
supported basis functio ns is achieved in the same manner 

although the book-keeping is 



as described in Section 2.2 



more complicated. Similarly the partial Householder trans- 
form required to separate the / = and 1 = 1 modes need 
only be applied to the m = and m = ±1 blocks of the 
resultant conversion matrix (Appendix An important 
side-effect of the separation in m is that the Yl(r) have the 
same trigonometric ^-dependence as the full-sky spherical 
harmonics (Appendix ^). This also implies that the Yl{r) 
can be treated as two-index quantities, defined by m and a 
second, arbitrary index in place of I. 

The algorithms described here were implemented on the 
Cambridge Centre for Mathematical Science's COSMOS 64- 
processor Silicon Graphics Origin 2000 and the evaluation 
and decomposition of the coupling matrix at the highest 
Planck resolution of imax — 2500 required about an hour. 
The majority of the time was spent factorizing the sub- 
matrices, and thus significant accelerations are unlikely, the 
highly-optimised linear algebra package (lapack; An- 
derson 1992) routines having been used for all the decompo- 
sitions. For a given choice of 6\ and 62, the decomposition 
of C need only be performed once, so orthogonalisation of 
the spherical harmonics on an incomplete sky should com- 
prise only a small fraction of the analysis required for the 
forthcoming MAP and Planck missions. 



3 HARMONIC ANALYSIS 

Methods for constructing an orthonormal basis set on the 
incomplete sphere from the spherical harmonics were dis- 
cussed in Section ^ but in most cases of data analysis it 
is the harmonic coefficients, representing functions on the 
sphere, that are of interest. There are at least three useful 
harmonic expansions of a general function on the sphere, 
and the relationships between these coefficients, which are 
summarised in Table ^, are derived here. 

A band-limited function, a(f), can be completely spec- 
ified by a finite number of harmonic coefficients as (cf. 
Appendix ^) 



a{r) = Y (r)a, 



(19) 



where it is assumed that imax is greater than or equal to the 
band-limit of a(r) and the harmonic coefficients are defined 
by 

a= Y(r)a(r)dn. (20) 
Js 

The invertibility of these transformations is due to the or- 
thonormality of the spherical harmonics on S and the fact 
that they represent a complete basis set given the band- 
limit. 

If a(f ) is only known over some fraction of the sphere 
S' < S, then a cannot be determined as above, as the inte- 
gral in equation ( pc| ) is incomplete. In this case the psuedo- 
harmonics 



Y{r)a{r)dn 



(21) 



fully specify a{r) in S' due to the band-limit. From equa- 
tion (Um they are related to the full harmonic coefficients 
by 



Ca, 



(22) 



where C is the coupling matrix, defined in equation (|l|). 

The psuedo-harmonics are useful quantities, but it is 
preferable to work with basis functions that are orthonor- 
mal on S' . Denoted Y'{r) in Section ^, their harmonic co- 
efficients are given by 



Y'{r)a{r)dn. 



(23) 



The relationship between Y{r) and Y'{r) given in equa- 
tion (^ fiows through to the harmonic coefficients and ap- 
plying equations (^), (|l^) and (^9|) to equation ( p^ ) gives 



's 

= BCa 



B / Y{r) [Y'^{r)a] dQ 
Js' 



(24) 



2.2) 



but, 
the non- 



The form of the above transformation depe nds on the 
decomposition used to generate A (cf. Section 
for CMB analysis, it is desirable to separate 
cosmological modes. This amounts to demanding that the 
only four of the have any contribution from the I = and 
I = 1 spherical harmonic coefficients. The uppper triangular 
structure of A""" = L as generated by a Cholesky decomposi- 
tion inherently satisfies this requirement, but in general the 
conversion matrix needs to be transformed explicitly. One 
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Table 1. Conversions between harmonic coefficients 



a = y(f)a{f) dn a=BTBa^C-ia a = B^a' ^ (A^) ^ a' 

d= \ Y{r)a(r)dQ, a = Ca d = Aa' 

a' = /g, Y'{r)a{r) dU a' = A^a a' = Bd 



The conversions between the various harmonic coefficients defined in Section ^: a are the 
standard coefficients o '_t le spherical harmonics [equation ( ^o| )] ; d are the psuedo- larmonic 
coefficients [equation (21)]; a' are the cut-sphere harmonic coefficients [equation (2^]; and 
d are the reconstructed spherical harmonic coefficients. The two basis functions are the 
spherical harmonics, Y{r) [defined in Appendix and the orthogonalised harmonics, 
y'(f) [defined in Section eI. If the coupling matrix of the spherical harmonics, C, is invert- 
ible, then the expressions for d following the can be used as exact inversions; otherwise 
the 'estimators' are approximate projections onto the cut region. A can be any matrix such 
that AA^ = C, and B can be any matrix such that BA = I. 



option is to use successive partial Householder transforms, 
as described in detail in Appendix ^. As can be seen from 
Fig. ^ the index-ordering and decomposition method com- 
bine to give a wide variety of conversion matrices; which of 
these is most suitable depends on the application. 

It is also possible to convert between the psuedo- 
harmonics and the cut-sky harmonics, as they both con- 
tain information about a{r) in 5" alone. Combining equa- 
tions (p^), ( pl| ) and ( [23I ) implies that a = Aa' and a' = Ba. 
However it is not always possible to determine a from either 
a' or a. In the low-Zmax limit these inversions are defined 



(Section 3.1), but for appreciable band-limit s on ly a projec- 
tion onto the cut sphere is possible (Section 3.2). 



3.1 Low-resolution analysis 

If the coupling matrix is numerically non-singular (i.e. 
^max ^ 50) then equation (bj) can be inverted to give 



and equation (ji^) implies that 
a — L a. 



(25) 



(26) 



These are specific examples of the fact that a band-limited 
function is completely defined if it is known over any fi- 
nite portion of the sphere, and a cut-sky analysis serves no 
purpose - any apparently localised contaminants infect the 
entire sky. However any measurement of a field on the sphere 
is subject to noise which is not band-limited, in which case 
the application of a cut has the effect of greatly amplifying 
the noise in the removed region (see Section 4.4), justifying 
the use of a cut-sky analysis in the low-resolution limit. 

3.2 High-resolution analysis 

If 'max <^ 50 and the coupling matrix is numerically sin- 
gular, it is impossible to reconstruct (even) a band-limited 
function that is known only on S' . The loss of information 
about modes constrained to the cut makes it clear that the 
analysis has the desired effect of removing contaminated (or 
otherwise problematic) regions, but the most appropriate 



transformation from the cut-sphere basis to conventional 
harmonics is less obvious. 

A least squares-like approach leads to a definition of the 
reconstructed full sky coefficients as 



(27) 



Similarly, equation (|2^) implies that 

d ~ B'^Ba = (AB)'^a, (28) 

where (AB)""^ is a projection operatorp] onto the range of 
C, which in real space is S' . If 'max —> 00 it is possible to 
write a{r) — ■ws'{r)a{r), where wg/(f) is the sharp window 
function defined in equation (^). 

It is at this point that the subtle distinctions between 
the use of a discrete cut and a window function become ap- 
parent. These results only hold for band-limited functions, 
but, as defined above, a{r) is not band-limited, and so can- 
not be analysed self-consistently. Whether a discrete cut or 
an apodizing function is to be preferred depends on the sit- 
uation in which the incomplete sky analysis is required. 



4 CMB DATA ANALYSIS 

In order to determine the properties of the C MB from noisy 
observations of the microwave sky (Section [4.l[ ) a number 
of non-trivial analysis steps are required, including: map- 
making (Section [4.2| ); component separation (Section 4.3) 
and power spectrum estimation (Section [4.41 ). Several algo- 
rithms have been suggested for all these steps, and these 
are discussed briefiy below, but the main focus is on when 
and how to apply a sky cut. Further, whilst it is possible to 
analyse the data in either real or Fourier space, the latter 
approach is emphasised here as it is more directly related to 
the formalism described in Sections ^ and ^ as well as being 
the focus of a related series of papers (van Leeuwen et al. 
2001; Challinor et al. 2001; Stolyarov et al. 2001). Note that 
the term 'map' is used here to denote any representation of 



^ The definition of B given in Section 2.2 implies that [(AB)'^]^ : 
(AB)"'^ and it is hence a projection operator if I 

max ^ 00. 
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i{l,m) 



i{l,m) 



i{l,m) 






i{l,m) 



i{l,m) 
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i{l,m) 







(h) 












X. 








i{l,m) 



Figure 3. Several examples of the conversion matrices, A^, that relate the orthonormal cut sphere harmonic coefficients, a' to the 
conventional spherical harmonic coefficients, a, according to equation (^). In all cases a symmetric, constant latitude of 6cut = ±20 
deg has been applied and imax = 10. The row index, i' , corresponds to the cut-sky harmonic coefficients and the column index, i{l,m), 
corresponds to the original spherical harmonic coefficients. The colour-map shows the absolute value of the elements of varying from 
zero (white) to their maximum value (black), which is normalised separately in each case. The spherical harmonic coefficients are indexed 
using i-ordering in the left column [(a), (d) and (g)]; m-ordering has been used in order to utilise the decoupling of different m-modes 
resulting from the azimuthal symmetry by using block-by-block storage in the central column [(b), (e) and (h)]; and m-ordering is used, 
but the structure of the coupling matrix is not utilised, in the right column [(c) and (f)]. In the case of (-ordering the non-cosmological 
modes are the four left-most columns of the conversion matrix; in the case of m-ordering the non-cosmological modes are the left-most 
columns of the three central blocks. The conversion matrices shown in the top row [(a), (b) and (c)] result from a Cholesky decomposition; 
those in the middle row [(d), (e) and (f)] are produced by an SVD (in which all modes with Wi^i < 0.1 have been removed); finally, in 
the bottom row [(g) and (h)], the partial Householder transform described in Appendix |c is applied to the conversion matrices shown 
in the above panels [(d) and (e), respectively]. 



a field on the sky and can imply either a set of real space 
pixel values or a vector of spherical harmonic coefficients. 



4.1 Observations 

Observations of the CMB can be made using a number of 
quite distinct techniques. Data have been obtained from the 
ground, high altitude balloons, and satellites, but the more 
important distinction is the type of telescope. The exper- 
iments listed in Section ^ include: straightforward single 
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dish telescopes, such as BOOMERanG (e.g. Netterfield et 
al. 2001) and Planck (e.g. BersaneUi et aL 1996); differenc- 
ing experiments, such as COBE (e.g. Smoot et al. 1992) 
and MAP (e.g. Jarosik et al. 1998); and interferometers, 
such as the Cambridge Anisotropy Telescope (CAT; Scott 
et al. 1996), the Cosmic Backround Imager (CBI; Padin 
et al. 2001) and the Degree Angular Scale Interferometer 
(DASI; Halverson et al. 2001). The interferometry surveys 
inevitably cover only a small fraction of the sky, and so a flat- 
sky Fourier analysis becomes possible. However the both the 
differencing and single dish surveys can, in principle, cover 
most of the celestial sphere, and should yield maps of the 
microwave sky that are limited only by (the combined effects 
of) instrument noise and the finite telescope beam. 



4.1.1 Noise 

The typical receivers used in the above experiments have 
two main noise contributions: random white noise and a cor- 
related low-frequency (i.e. '1//') component. The latter is 
potentially troublesome, leading to 'stripy' maps with cor- 
related errors, and is the main reason for the popularity of 
differencing experiments which remove low frequency noise 
at the moment of observation. However data from single dish 
surveys can be 'de-striped' if the scan strategy includes suf- 
ficiently many multiply-observed points (e.g. Tegmark 1997; 
Delabrouille 1998; Maino et al. 1999) or the time-time noise 
covariance matrix can be fully included in the map-making 
process (Wright, Hinshaw & Bennett 1996; Natoli et al. 
2001; Challinor et al. 2001). Hence correlated errors are ig- 
nored in the simple analysis presented here. 

This leaves only the white component, which can be 
analysed most simply in the case of a single beam exper- 
iment. Following Knox (1995), a receiver is characterised 
by its sensitivity, s (generally chosen to have units of tem- 
perature time^^^). Assuming the noise is Gaussian it has 
expectation values (n) — and (n^) — s^t over an integra- 
tion time t. The manner in which this noise projects onto 
a sky map depends on the map-making algorithm, the scan 
strategy, and the beam. 



4.1.2 Beam convolution 

All telescopes necessarily have a finite point-spread function 
or beam, which, for a given detector, can be characterised 
by b{r), the fraction of photons from direction r that are 
registered, given a nominal orientation towards the north 
pole (i.e. 9 = 0). The harmonic expansion of the beam in 



m 



this orientation is denoted b = hi 



{l,m) 



with the band-limit 



being related to the nominal resolution of the detector. For 
a given type of telescope the resolution improves with fre- 
quency due to diffraction effects; this places limitations on 
the component separati on a lgorithms that are used on the 
incomplete sky (Section 4.3). 

Most experiments have beams that are manifestly 
asymmetric, a fact which must be accounted for explicitly 
by the data analysis algorithms, but the cut-sky issues of 
interest here can be explored more clearly if the beam is 
approximated by its azimuthally averaged counterpart (e.g. 
Challinor et al. 2001). Defined by 



27r 



(29) 



its harmonic coefficients are simply 

bl.m = Sm.obl.O 

= 5^,0 y/ [21 + 1)71 [ b [arccos(x)] Pi {x)dx, (30) 



where Pi{x) is a Legendre polynomial (Appendix]^). The 
use of b{r) allows the definition of a beam-smoothed sky, 
s(f), given in terms of the true sky, s(f), by 



s(f) 



J b [arccos (r ■ f ')] s(f ') dSl'. 



(31) 



This convolution i s mu ch simpler in harmonic space, and 
applying equation (A6) to equation (^) yields 

s = Bs, (32) 

where B (as distinct from the conversion matrix, B) is a 
diagonal 'convolution matrix' with 



47r 



2l(i) 



Ol(i),0 



= 27r 



6 [arccos(a;)] P;(i) (x) dx, 



(33) 



l{i) being defined in Appendix^. The simple form of equa- 
tion (|3^) is often utilised explicitly in CMB analysis algo- 
rithms (e.g. Knox 1995; Hobson et al. 1998; Oh et al. 1999; 
Stolyarov 2001), but in all cases full sky coverage is - some- 
times implicitly - assumed. 

Turning to convolution on the incomplete sphere, S' , 
application of equation (^|) to equation (^^ yields 

s = A^Bs, (34) 



where A is defined implicith 
resolution limit equation (Gf 
logue of equation ( p^ as 

s = B s , 



in equation (^^. In the low- 
) then gives the cut-sky ana- 

(35) 



with the new convolution matrix defined by B' = A^'^BB''". 

For higher band-limits no such relation exists, the loss 
of modes in the cut region rendering the convolution ill- 
defined. This is simply understood in real space, as the value 
of s(r) near the edge of S" is given by an integral that ex- 
tends several beam widths into the removed region. Thus it 
is impossible to relate s(f) to s(f) with r G S' . These ar- 
guments are true independent of the representation chosen, 
but in harmonic space they mean that it is impossible to 
relate s' to s' . 

Whilst equation ( psj ) is formally incorrect in high-Zmax 
cases, it is potentially useful as a practical approximation. 
It is equivalent to assuming that the signal is given by 
Y'{r)s' , which implies that s{r) ~ in S — 5". This is par- 
ticularly inaccurate if the removed region contains anoma- 
lously strong sources, such as the Galactic plane. Nonethe- 
less, equation (|35) gives s'(f) correctly for all r more than 
a few beam widths away from the edge of the cut region. 
However, even if this is an acceptable approximation, there 
is the further inconvenience that the effective cut-sky beam, 
B , is not diagonal, introducing couplings between all the 
modes. 
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In short, it is preferable to avoid performing any sort of 
convolution (or deconvolution) on the cut sky, although it is 
clear that this situation is encountered in any survey with 
incomplete sky coverage. The one, albeit trivial, exception 
to this rule is if the beam is a delta function, or at least 
the closest approximation to a delta function possible given 
the band-limit under consideration. In this case B = I and 
hence, from equation (p^), B' = T as well. Equation ( |35[ ) 
then implies that s' (= s') is the true, unconvolved sky 
map, estimation of which is addressed next. 



4.2 Map-making 

Some of the most important products of the next genera- 
tion of CMB survey will be high-resolution maps of the sky 
at each of several frequencies. Such maps can be created in 
a number of ways, but care must be taken to account for 
a huge variety of systematics whilst retaining as much in- 
formation as possible. Both real space (e.g. Wright et al. 
1996; Bond et al. 1998; Natoh et al. 2001) and Fourier space 
(van Leeuwen et al. 2001; Challinor et al. 2001) algorithms 
have been proposed as being suited to particular apsects 
of the map-making problem. The resultant uncertainties in 
the pixel values or harmonic coefficients depend on both 
the data itself (i.e. the scan strategy, noise properties, etc.) 
and the map-making algorithms used, and can vary quite 
markedly from experiment to experiment. 

Here only the idealised case of uniform sky coverage is 
considered as the discussion which follows is not significantly 
changed by this useful simplification. Under this assumption 
the optimal estimator for the unsmoothed sky, s {= s -\-ri), 
would be unbiased (i.e. (s) = s) and have covariance given 
by 



N 



{n'n?'^ = <^(s — s) (s — s) 



1 



a a , 



(36) 



where (cf. Knox 1995) = 47rs^/(A'dfobs), iobs is the to- 
tal observation time of the survey, and A'^d is the number of 
detectors at the frequency in question (all of which are as- 
sumed to have the same beam). An important issue at this 
point is the band-limit chosen. Clearly only a finite analysis 
is possible in practice, and B becomes increasingly singular 
as /max oo; these two points are related in so far as the 
sky can never be reconstructed with infinite resolution. The 
choice of imax is somewhat arbitrary, although any value a 
factor of a few greater than the effective beam width will en- 
sure that B exists whilst discarding only multipoles that 
are noise-dominated. The fact that, unlike the useful sig- 
nal, the noise is not subject to any band-limit is critical to 
the understanding of the low-^ max cut-sky power spectrum 



estimation discussed in Section 4.4. 

Note also that, due to the assumption of uniform sky 
coverage, the covariance matrix is diagonal. Transforming 
this estimator into real space yields maps with covariance 
given by 



{n{ri)n{r2)) 



(37) 



As the noise term in the data is not beam-convolved the 
removal of the beam results in spatial correlations of the 
noise (as encoded in B), as well as correlations due to the fi- 
nite resolution analysis (the sums over spherical harmonics) , 
which are essentially equivalent to pixel smoothing. In the 



more realistic case of non-uniform sky coverage, the covari- 
ance matrix is non-diagonal in both bases, a point discussed 
further by Oh et al. (1999). 

The above estimator for the true sky is closely linked to 
the more commonly used estimator for the smoothed sky, s. 
Being related by s = B s, it is clear they contain the same 
information (under the assumption of a symmetric beam). 
The covariance structure of s is simpler as the correlations 
discussed above are not introduced, but s is a more natural 
data object in the context of this discussion as it is the true 
sky that is of interest. In particular, unsmoothed maps al- 
low more flexibility in applying a sky cut, as the problems 
with convolution on the incomplete sphere described in Sec- 
tion 4.1./ do not arise. In practice the best compromise may 
be to reconstruct the sky convolved with the azimuthally 
averaged beam, thus creating maps with the simplest co- 
variance structure possible without information loss. This 
can be done in either real space (e.g. Bond et al. 1998) or 
harmonic space (Challinor et al. 2001), although the real 
space pointing matrix is more complicated if beam asym- 
metry information is included. 

In summary, if a survey covers the entire celestial sphere 
it is preferable to use full-sky frequency maps. However it 
is possible that small parts of the sky will be missed due 
to either the scan strategy (cf. Maino et al. 1999) or hard- 
ware problems during the survey itself. If this is the case 
the best unsmoothed map that could be constructed would 
be larger than the actual observed region, but the errors 
around the boundary of this area would be very high. An 
inferential approach is possible, but significant difficulties 
are encountered, especially in Fourier space (Challinor et al. 
2001). Fortunately, it is probable that both MAP and Planck 
will produce full sky maps at several frequencies, which can 
then be used to construct maps of the various astrophysical 
components. 

4.3 Component separation 

The microwave sky consists of several distinct astrophysical 
components, as listed in Section |l[ Fortunately they have 
sufficiently distinct spectra that they can be separated us- 
ing multi-frequency data. Given that MAP and Planck will 
produce maps in five and ten bands, respectively, it should 
be possible to produce maps of the various components (par- 
ticularly the CMB) that are relatively free of contamination. 
As with map-making, a number of algorithms have been put 
forward for this stage of the data processing, although the 
main focus has been on Fourier space methods (e.g. Tegmark 
& Efstathiou 1996; Hobson et al. 1998; Bouchet & Gispert 
1999; Prunet et al. 2001; Stolyarov et al. 2001). Aside from 
the expected statistical isotropy of the CMB signal, one of 
the reasons for this emphasis has been th e simp licity of beam 
convolution in harmonic space (Section 4.1.2). This is crit- 
ical if smoothed frequency maps are used as the effective 
smoothing scale will vary with frequency if the telescope is 
(close to) diffraction limited. However if unsmoothed maps 
are used real space component separation methods (e.g. Bac- 
cigalupi et al. 2000) must also come into consideration, the 
optimal choice of basis being less clear. 

One common aspect of all the separation techniques 
referenced above is that much of the (prior) information 
about both signal and noise correlations is disregarded in or- 
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der to render the problem computationally feasible. In real 
space the correlations between nearby pixels are ignored, 
and in Fourier space it is the mode-mode couplings that 
are neglected. Surprisingly, these approximations appear to 
be unimportant in practice - even the Galactic components 
have been recovered with striking accuracy. The most rele- 
vant result to this discussion is the all-sky component sepa- 
ration to Zmax — 2500 presented by Stolyarov et al. (2001), 
as it provides clear evidence that whatever correlations are 
present in the full-sky harmonic basis are unimportant - 
there are some errors close to the Galactic centre, but they 
are localised, and there is no sign of this affecting the recon- 
struction globally. 

If the Galactic plane is removed prior to component 
separation this one troublesome region is no longer present 
in the analysis, but new problems arise. Firstly, smoothed 
maps (with frequency-dependent beam-widths) cannot be 
used as input data without inducing errors around the edges 
of S" due to the ill-defined nature of convolution on the cut 



sky (Section 4.1.2). Even if such errors are deemed accept- 
able (e.g. Prunet et al. 2001) or beam-deconvolved maps are 
used, the transformation described in equation (^^ com- 
pletely changes the correlation structure of the harmonics. 
In particular, the signal-signal correlation matrices are non- 
diagonal for all components, including random fields like 
the CMB (see Section 4.4). Whereas the couplings between 
the spherical harmonic coefficients can apparently be disre- 
garded, this has not been demonstrated for these induced 
correlations in the orthonormal basis. Prunet et al. (2001) 
performed cut-sky component separation including them in 
full, but were thus limited to Zmax — 500, the computational 
task being made considerably larger. 

In real space the application of a cut is trivial, pro- 
vided that beam-deconvolved maps are used, as it simply 
requires that pixels in the removed region be ignored. Thus 
the Baccigalupi et al. (2000) method should be well suited 
to a cut-sky analysis. 

Given that realistic component separation simulations 
have only recently become available, it is likely that im- 
portant developments in this field will be made in the near 
future. For the moment, however, it appears that separa- 
tion can usefully be performed on either the full or cut sky 
without introducing catastrophic errors. Provided the mod- 
els of microwave emissions from the Galactic plane used in 
the above simulations are sufficiently realistic, it may thus 
be preferable to generate the full-sky maps of the various 
astrophysical components, retaining the option of masking 
unwanted regions at a later stage. 

4.4 Power spectrum estimation 

If the fluctuations in the early universe were the result of 
inflation (e.g. Linde 1990) then the CMB is expected to be 
a Gaussian random field, the statistical properties of which 
can be specifled completely by its angular power spectrum, 
Ci. Even if this is not the case, the power spectrum should 
encode much of the cosmological information present. It is 
thus unsurprising th at, a s with map- making and component 
separation (Sections 4.2 and 4. J, respectively), many differ- 



ent methods of power spectrum estimation have been de- 
veloped (e.g. Tegmark 1997; Gorski 1994; Bond et al. 1998; 
Oh et al. 1999; Szapudi et al. 2001; Wandelt et al. 2001; 



Hivon et al. 2001). Further, sky cuts have been incorpo- 
rated into many of these algorithms as it seems certain that 
the strength of the Galactic microwave emissions will pre- 
vent the CMB from ever being accurately measured in this 
region. Due to the proliferation of papers on this subject, 
this discussion of power spectrum estimation is limited to 
a description of a maximum likelihood formalism using the 
orthonormal basis functions described in Section ^, with ref- 
erence to how their behaviour differs in the low- and high- 
resolution regimes. 



4 ■4-1 Maximum likelihood formalism 

The most powerful method of power spectrum estimation is 
maximum likelihood (e.g. Press et al. 1992), although this 
has only been implemented to MAP resolution to date (Oh 
et al. 1999). By invoking Gaussian statistics for both the 
CMB and the noise, it is possible to write down the exact 
likelihood for the observed map (e.g. Gorski 1994; Borrill 
1999). On the full sky the effective data- vector is s = s + n 
(i.e. the est imator for the true sky, of the form discussed in 
Section 4.2 and not the quantity being estimated here) with 
the non-cosmological I = and I = 1 modes removed. The 



full likelihood is given by 

dp _ 1 
dl ^ 



7(2^ 



det(S + N) 



: exp 



is^(S-f N)-is 



(38) 



where S and N are the signal and noise covariance matri- 
ces, respectively. The assumption of Gaussianity implies that 
•5*;,^' = SiyCu^i), where C; is the CMB power spectrum and 
l{i) is defined in Appendix The form of N is determined by 
a combination of the survey method and the data processing 
up to this point, but is unlikely to have the simple form of 
equation (^) due to the imperfect component separation. 
The maximum likelihood calculation consists of finding an 
estimator for the underlying power spectrum, Ci , such that 
equation (^) is maximised, and there are a number of algo- 
rithms for finding this quantity (e.g. Bond et al. 1998; Oh 
et al. 1999). 

The maximum likelihood formalism on the cut sky takes 
the same form as on the full sky, but with the data-vector 
s' = A^s and the covariance matrices suitably transformed 
to give (cf. Gorski 1994) 

dp _ 1 



exp 



max-4det(S' + N') 

-s (S + N ) s 



(39) 



with the four modes con taining information on the monopole 
and dipole (Section 2.2) again excluded. The signal covari- 
ance matrix is subject to a simple similarity transform, 
S' = A^SA, but the same is not true for the noise covariance 
matrix as the noise field is not band-limited (a fact critical 
to the use use of a cut-sky analysis at low resolution, as 
discussed below) . The coupling of the cut-sky modes makes 
the maximisation of equation ( p9[ ) non-trivial (cf. Oh et al. 
1999), even if S -I- N is diagonal on the full sky. Nonetheless 
it is useful to work under this idealised assumption in or- 
der to see how the application of the cut ensures that the 
maximum likelihood solution is independent of the Galactic 
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signal; the manner in which this is achieved is quite difTerent 
in the low- and high-resolution cases. 

^.4-2 Low-resolution analysis 

The effect of a sky cut on power spectrum estimation is not 
entirely obvious in the low-Zmax case in which the coupling 
matrix (Section 2.1) is invertible. The effective band-limit, 
produced by th e co mbined effects of the beam and noise 
(Sections ^ and 4.2), means the signal over the whole sky 
(including e.g. the Galactic plane) is encoded in the cut-sky 
coefficients. Thus the application of a cut would be redun- 
dant were it not for the presence of non-band-limited noise 
which cannot be characterised properly a finite harmonic 
analysis. 

Applying the incomplete spherical transform defined in 
equation ( ^jj) to a purely white noise field n(f) [i.e. (n(r)) — 
and {n{ri)n{r2)} = 5{ri — f2)(J^; cf. equation ([37])] gives 
cut-sky harmonic coefficients with (n') — and 

N' = (n'n'^) = o-'l. (40) 

Projecting back into real space gives a field n'(f) which sat- 
isfies (n'(f)) = and 

{n'(fi)n'(f2)) =a"y'^(fi)C-^y(r2). (41) 

On the full sky the same procedure (i.e. a finite spherical 
harmonic analysis followed by a transformation back into 
real space) would yield a noise field with covariance structure 
given by 

{n(fi)n(f2)) = (7"Y^(fi)Y(f2) 

2 

= ^T.^^l + miri-r2), (42) 

where Pi{x) is a Legendre polynomial (Appendix^. Taking 
the limit r2 — > ri this implies that {n^{f)) — imaxO"^/(47r), 
which represents smoothing relative to the original noise 
field caused by the use of a finite analysis. 

Whilst this smoothing occurs in both the cut- and full- 
sky formalisms, the presence of in the former case [equa- 
tion (Ull)] implies a spatial dependence. As can be seen from 
Fig. W^he noise in the cut region is greatly increased, which 
is a natural way of formally encoding the qualitative fact 
that, for whatever reason, the data in the cut is contam- 
inated by more than just the white noise field. Thus, de- 
spite the invertibility of the coupling matrix (and the band- 
limited cut-sky analysis), the application of a cut has the 
desired effect of greatly reducing the impact of any spurious 
signal, such as the Galaxy. However Fig. ^ also implies that 
a similar effect could be achieved without performing a cut 
(and hence leaving the signal unchanged), instead adding a 
high level of artificial noise in the offending region(s). Fi- 
nally, it is important to note that the dual assumptions 
used in the derivation of equation ( pl| ) - uniform noise and 
no beam - are unrealistic, but the manner in which a low- 
resolution cut-sky analysis works is the same in less idealised 
scenarios. 



4-.4-3 High-resolution analysis 

The high-resolution case is more straightforward, as the ap- 
plication of the cut results in a data- vector, s', which con- 



tains little information about the removed region. This is 
quite distinct from the low-resolution case discussed above, 
in that here it is the predominantly the signal that is 
changed, rather than the noise. That said, the noise close 
to the boundary of the cut is increased in the same manner 
as explained above. This has the same effect as the apodiz- 
ing function formalism described by Tegmark (1997), down- 
weighting points around which there is not full correlation 
information. 

Another difference between the low- and high- resolu tion 
analyses is that s' is smaller than s, from Section [2.1.l| . Al- 
though this does not result in any significant computational 
saving, it serves to emphasise the information loss associ- 
ated with removing part of the sky, and is an independent 
derivation of the fact that the uncertainties in the estimated 
power spectrum increase as Att/^Is' (cf. Hobson & Magueijo 
1996; Tegmark 1997). 



5 CONCLUSIONS 

The upcoming microwave surveys will require a cut-sky 
analysis to prevent the strong Galactic emissions from con- 
taminating the CMB signal. The spherical harmonics are 
non-orthogonal on the cut-sphere, but an orthonormal basis 
set can be constructed from them using SVD-based tech- 
niques (Section 1^). The application of the resultant conver- 
sion matrix to the conventional multipoles results in cut- 
sphere harmonics that contain only the desired information. 
In the low-resolution case the influence of the Galaxy is re- 
duced by increasing the effective noise in the cut; in the 
high-resolution limit the orthonormal basis functions can 
model the infinitely sharp cut sufficiently well that they 
have no support in the removed region. It is also important 
to note that the cut should probably only be applied after 
beam-deconvolution has been attempted, as convolution is 
ill-defined on the incomplete sphere. 

The algorithms described here were implemented to 
Legendre multipoles of /max — 2500 for a constant lati- 
tude cut, in which case the coupling matrix of the spher- 
ical harmonics is block-diagonal. At present computational 
limitations make a general orthogonalisation impractical for 
/max 200, although there are some possibilities to extend 
this. For instance only ~ 1 per cent of the coupling ma- 
trix contains significant information if the cut is well-chosen 
(e.g. rectangular in 9 and </)) and so sparse matrix techniques 
should thus allow orthogonalisation to /max — 1000 in this 
case. 

Another requirement is orthogonalisation of tensor ba- 
sis functions on the incomplete sphere, as both the MAP and 
Planck satellites will measure polarization. The resultant 
formalism is more complicated, but the same general princi- 
ples hold; this issue is explored further by Lewis, Challinor 

6 Turok (2001). 
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Figure 4. The variance of a ("in itia lly) uniform Gaussian noise field as a function of latitude, b, after application of the cut-sky orthog- 
onalisation described in Section Constant latitude cuts of bcut = 10 deg (a) and bcut = 20 deg (b) were applied (as indicated by the 
dashed vertical lines) and results are shown for /max = 20 and (max = 40, as labelled. (The oscillations near the peak of the latter curve 
are indicative of the limited accuracy of the decomposition of the ill-conditioned couling matrix.) 
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APPENDIX A: SPHERICAL HARMONICS 

The spherical harmonics form a complete set of orthonormal 
basis functions over the entire sphere. They are most com- 
monly defined as complex functions (e.g. Landau & Lifshitz 
1976; Brink & Satchler 1993), but it is more convenient to 
use real harmonics in this application. Adapting the nota- 
tion of Gorski (1994), the real spherical harmonics are given 
by 



Yi,^(r) = Yi,™(e',0) = A,j,„|[cos(6»)]s,„(<^), 
where I > and |m| < I and 

%/2 sin(|m|<jf>), if m < 0, 

Sm(<^) = { I, if m = 0, 

\/2 cos(|m](j!)), if m > 0, 



(Al) 



(A2) 



implying that " Sm(<^)s^/ (0) d(j} = 27r(5„_„/ . For < m < 
I and — 1 < a; < 1 the normalised associated Legendre func- 
tions are defined by 



2l + l(l-m)\ „ , ^ , . , 

A,.(.) = ^^i^P,.(.). (A3) 

Hence J^^ A!,m(s)A;',m(K) da; = (5j ;//(27r). Under the same 
conditions the (unnormalised) associated Legendre functions 
are given by[J| 

p,™(.) = (-i)-(i-x^)'"/^iL_p,(,), 

with the Legendre polynomials given by 



(A4) 



(A5) 



2'1\ dx' 

A real field on the sphere, a{r), can be expanded in 
terms of spherical harmonic coefficients, given by 



II This definition, with the (— l)"" term, correpsonds to that 
given by Abramowitz & Stegun (1971) and Gradshteyn & Ryzhik 
(2000) but differs from that used by Arfken (1985) and Brink & 
Satchler (1993). 



fli.m = / Yi,m(r')a(f ) dn. 
Js 

This can be inverted to give 

!max I 
i = m=-i 



(A6) 



(A7) 



provided that Zmax oo, due to the orthonormality of the 
spherical harmonics on the full sphere: 



(A8) 



If a finite /max is used this inversion is no longer possible 
for general a(f), although it does still hold for band- limited 
functions^. 

Whilst the two indices I and m have quite distinct in- 
terpretations it is convenient to combine them into a sin- 
gle index, i, which allows the definition of vectors Y{r) = 
Yi(i,m){'f') and a — ani^^y Two obvious indexing schemes 
present themselves: grouping in / and m. The first, as in- 
troduced by Gorski (1994), is natural for power spectrum 
estimation and very simple: 



i{l,m) = I +l + m+l. 

The two 'inverses' of this relationship are 

/ = int [(i-1)'''''] 

and 

m ^ i~ {[■^ + I + 1). 



(A9) 
(AlO) 
(All) 



The second choice of ordering is useful in cases of azimuthal 
symmetry in which the orthogonality expressed in equa- 
tion (A2) is maintained, and grouping in m is achieved by 
defining 



m) = 



l + m + 1 

+ (imax + m)(/max + m -I- l)/2, if m < 0, 



max max + 1)^ 

- (^max - m)(/max - m + l)/2, if 771 > 0. 

(A12) 

The 'inverses' in this case are given by 

'int [-«n,ax + i (\/8iTT - 3)] , 

if l< (;max + l)(/max + 2)/2, 



mt 



max + 1)2 - I + 1] + 1 - 3j 
if i > (?max + l)(/max + 2)/2 

(A13) 



and 



(A14) 



i — [m + 1 

+ (^max + m)(;nmx + m + l)/2] , if m < 0, 
i - [-/max + {I max + 1)^ 

-(/max - m)(Zmax - m + l)/2] , if 771 > 0. 



** A band-limited function can, by (somewhat circular) defini- 
tion, be constructed from a finite sum over spherical harmonics. 
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Other indexing schemes have been used in the more specific 
case of simuiated Planck data-sets in which the slcy cover- 
age is periodic in azimuth (van Leeuwen, private communi- 
cation), but are beyond the scope of this paper. 



APPENDIX B: INTEGRATION OF PRODUCTS 
OF ASSOCIATED LEGENDRE FUNCTIONS 



In Section 2.3 



integrals of the form 

\,m{x)\v ^rn{^) da; 



(Bl) 



arose; here the \i,m{x) are the norm alis ed associated Leden- 
gre functions, defined in equation (A3), and m is assumed 
to be non-negative. These integrals can be evaluated quickly 
and accurately using a combination of closed formula and 
recursion relations. 

The a ssoc iated Legendre functions, Pi,m{x) [defined in 
equation (A4)], are solutions of the ordinary differential 
equation (e.g. Arfken 1985) 



_d_ 

da; 



2\ dPi,m 



+ 



Pi,m (a;) = 0.(B2) 



Multiplying this equation by P;/ „,(a;) and integrating (from 
xi to X2) by parts twice yields 



{l-l'){l + l' + 1) / Pi,,„{x)Pr^,{x)dx 

•J X\ 

(1 - x') p^^{x)^^ - (1 - x^) P,,^{x) 



(B3) 



dPi,„ 
dx 



This is a reflection of the standard result that integrals 
of solutions of a self-adjoint differential equation [as equa- 
tion ( B2 ) is] can be expressed as boundary terms (e.g. Arfken 
1985). The derivatives in equation (B3) can be removed by 



using the standard recursion relationship (e.g. Gradshteyn 
& Ryzhik 2000) 



(1 - X^) = + m)Pl-l,m{x) - lxPl,m{x) 

to yield, for I 7^ I' , 

^2 1 
Pi.m{x)Pr ^^{x) dx 



(B4) 



X [{l' +m)Pi,^{x)Pir_i^^{x) 

+ {l - l')xPl,m{x)Pi,,m{x) 
-{I + m)Pl.i^rn{x)Pl> ^rn{x)]\ 



X — X2 



(B5) 



Note that the first term must be omitted if = m and that 
the third term must be omitted \i I = m; these Legendre 
functions are implicitly zero from equation (|A4[) . Finally, 
this can be normalised according to equation (|A3l), giving 



2/ -f 1 
21-1 



{P - m2)Ai_i,™(x)Ai/_„(x) 



An alternative derivation of this result was presented by 
Wandelt et al. (2001); it is also in principle equivalent to 
Eq. 5.9(13) of Varshalovich, Moskalev & Khersonskii (1988), 
but their application of equation (B4) is in error. 

For the case I = I' a recursion relat ion is required, start- 
ing with I — m. Combining equations (A3) and (A4), 



\m,m{x) = (— 1)'' 



2m -I- 1 
4-K 



(2m- 1)!! (1 -a;^)' 



1/2 



(B7) 



where n!! = 1 x 3 x • ■ ■ x (71 — 2) x n for odd n. Integrating 
by parts and using equation (B7) again gives 



X2 —xi 
47r ■ 



'm — l,m — 1,771- 



-i(a;i,a;2) -I- 



(B8) 
if m = 0, 

if m > 0. 



Moving to Z = m -f 1, the standard relationship (e.g. 
Gradshteyn & Ryzhik 2000) that 



\m + l,m{x) — {2m + 3)xXm,m{x) 



combines with equation (B7) to give 
-''m+i,m-t-i,m(a;i, a;2) ~ 



ixi,X2) 

2m + 2 2 f x| 

^—^ X \^+,,^+,(x)\ 



(B9) 



(BIO) 



where the first term is given in equation (B8). 

The last step is to derive a recursion relation relating 

Il.l,ni{xi,X2) to /!_l,;_l,m(xi,a;2) and Il-2.l-2,m{xi,X2)- Eq. 

(C22) of Wandelt et al. (2001) gives a four-term recursion 
to obtain 7; ;/ „,(a;i, a;2); it can be applied successively (once 
swapping I and /') to obtain 



/;,i,m(a;i,a;2) = 



2m^ - 2r + 2? - 1 



(l -m){l + m) 
(;-l-m)(Z-l + m) , 



'!-l,i-l,m 



(a;i,a;2) 



{I -m){l + m) 

21-1 i 

2; + l {I - m){l + m) 



ii-2J-2,m 



(a;i,X2) 



^) h,m{x)\l-l^m{x)\'^^J^^ 



2l-l{l- 



m){l- 



21-3 {l-m)^{l + m)^ 

(1 - x^) A;_i,,„(a:)Ai_2,m(a:)r^'"^ ■ 

V / I x=xi 



In summary, equation ( pjq ) can be used to evaluate all 
Ii.i' . m(xi. X2) for which I ^ I , and equations (B8), (BIC) 



and (Bll) combine to give all /i,i,m(a;i, a;2) recursively. 



Il,l',m(xi,X2) 



{I - r){i + I' + 1) 



21' + 1 



{P - m2)Ai,,„(a;)A,/_i,„(a;) 



2/' - 1 ■ 

-I- (/ - l')x\l^m{x)\ii ^^{x) 



(B6) 



APPENDIX C: TREATMENT OF 
NON-COSMOLOGICAL MODES 

All the cosmological information encoded in the CMB is ex- 
pected to be contained in the I > 2 modes; the I — mode 
in an isotropic universe can be normalised arbitrarily and 
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the I = 1 modes can be set to zero by adopting an appro- 
priate reference frame. Nonetheless, observations of the mi- 
crowave sky will yield non-zero monopole and dipole values 
for a number of reasons (e.g. the observer's motion; Galac- 
tic emission; extra-Galactic point-sources). Hence these low- 
order modes must be included in the analysis of CMB data, 
but should be kept separate from the cosmological modes, 
as is naturally the case if spherical harmonic coefficients are 
used to describe the data. It is also important to note that 
the properties of the basis functions themselves are unim- 
portant - the essential requirement is that only four of the 
cut-sky harmonic coefficients contain information on the un- 
wanted modes. 

The method of orthogonalisation summarised in equa- 
tions (|ll|), and ( [l4| ) does not explicitly impose any 
particular structure on the conversion matrix, A [which re- 
lates harmonic coefBcients on the incomplete sphere to those 
on the full sphere by a' = A^a; equation (p4)]. The non- 
cosmological modes are kept separate from the cosmological 
modes if the first four columns of have only zeros from 
the fifth row on, assuming the full-sky harmonic coefficients 
are indexed using ^-ordering (Appendix^. This is achieved 
naturally if A^ is constructed to be upper triangular, as in 
the case o f the Cholesky decomposition described in Sec- 
tion |2.2.l[ The other decomposition methods discussed in 
Sectionpyd do not share this property, and so the resultant 
conversion matrices must be adjusted explicitly. 

One way of achieving this is to use a partial Householder 
transform (e.g. Press et al. 1992). The last i'^i^^ — i elements 
of the ith column of a general ij^ax x *max matrix can be set to 
zero by the transformation M' = PiM, with the orthogonal 
Householder matrix defined by 



applied to each of three blocks as is. The only slight incon- 
venience is that it is no longer the first four cut-sky modes 
that contain the non-cosmological information, and the rel- 
evant modes must be ffagged explicitly. 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science MJ^X style file. 



Pi = I 



(CI) 



where rrii is given by 



0, 



{mi)] = < 



M,- 



if j < i, 



if j > i, 



(C2) 



and i < mm{i'^^^, imax) is assumed. Provided that the 
Householder matrix applied to B is that generated from A""", 
the transformations A'""" — PiA^ and B' — P^B leave equa- 
tions ( p^ ) and ( p^ unaffected as Pi is orthogonal by con- 
struction. Applying Pi, P2, P3 and then P4 to the succes- 
sively updated A^ ensures that the I = and I = 1 modes 
influence only the first four cut-sky harmonic coefficients, 
as required. This procedure could be continued, moving A""^ 
successively closer to upper triangular form, although this 
cannot be achieved in full as A""" has more columns than 



Special mention must be made of the constant lati- 
tude cut case, the symmetry of which can only be utilised 
if the spherical harmonics are indexed using m-ordering 
(Appendix In this case only the m = and m = ±1 
blocks have any contribution from the monopole or dipole, 
and each can be treated separately. Further, the ordering 
within these blocks is such that the non-cosmological modes 
are in the first rows, and so the above algorithm can be 
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